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Three cellular automaton (CA) models of increasing complexity are intro- 
duced to model driven diffusive systems related to the generalized Frenkel- 
Kontorova (FK) models recently proposed by Braun et al. [Phys. Rev. E 

X" 

$— i ' 58, 1311 (1998)]. The models are defined in terms of parallel updating rules. 

Simulation results are presented for these models. The features are quali- 
tatively similar to those models defined previously in terms of sequentially 
updating rules. Essential features of the FK model such as phase transitions, 
jamming due to atoms in the immobile state, and hysteresis in the relation- 
ship between the fraction of atoms in the running state and the bias field are 
captured. Formulating in terms of parallel updating rules has the advantage 
that the models can be treated analytically by following the time evolution of 
the occupation on every site of the lattice. Results of this analytical approach 
are given for the two simpler models. The steady state properties are found 
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by studying the stable fixed points of a closed set of dynamical equations 
obtained within the approximation of retaining spatial correlations only upto 
two nearest neighboring sites. Results are found to be in good agreement with 
numerical data. 

PACS Nos.: 05.70.Ln, 05.45.-a, 66.30.-h, 05.60.-k 



2 



Typeset using REVTgX 



I. INTRODUCTION 



The physics of driven diffusive systems has attracted much attention recently [|l|,0 due 
to their relevance to the general area of nonequilibrium statistical mechanics and their wide 
range of possible potential applications. In particular, the Frenkel-Kontorova (FK) model 
HTJ] and its generalizations 0-0 have been studied within the context of tribophysics. Braun 
and coworkers |^-T^] studied the atomic current in one and two-dimensional atomic systems 
in the presence of a periodic potential under the influence of a dc driving force within 
the approach of Langevin equations. In tribophysics, the driving force emerges owing to 
the motion of one of the two substrates separarted by a thin atomic layer. The results of 
these studies are characterized by two features. One feature is that the system exhibits 
hysteresis in response to the driving force. The system jumps between low-mobility and 
high-mobility regimes in a hysteretic manner as a function of the driving force. Another 
feature is that accompanying this transition, the atoms tend to organize themselves into 
two types of domains consisting of atoms in states of different characters, one consisting of 
slowly moving ( "immobile" ) atoms and another consisting of "running" atoms moving with 



maximum speed. The latter feature resembles those in traffic flow models [TT|-|T3| in which 
cars may be moving at its maximum speed if they are not blocked or may be momentarily 
stationary if they are blocked in front. 

The models studied in Ref. |8HiO| are quite complicated. Attempts have been made |14 



to introduce simpler models which capture the essential features. In a series of three lattice 
gas (LG) models of increasing complexity (henceforth referred to as LG Models A, B, and 
C), Braun et al. introduced probabilistic hoppings of atoms on a lattice together with the 
possibility of the atoms being found in one of two possible states. The underlying model (LG 
Model A) in one dimension is that N atoms are placed in a lattice of length L corresponding 
to a concentration of p = N/L. The dynamics is introduced in a random and sequential 
fashion by randomly choosing a site at each time step and updating the system according 
to specific rules. The hopping to the nearest neighbouring sites of an immobile atom in a 
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randomly chosen site is characterized by a probability a (1 — a) that an atom hops into the 
site in the right (left) hand side. Thus the parameter a models the effect of a driving force, 
and the hoppings to the right and left hand sides correspond to the effects of drift together 
with diffusion. Atoms in the running state always attempt to hop to the right, which is 
taken to be the direction of the driving force. The a = 1/2 case corresponds to vanishing 
bias field. The a = 1 case corresponds to the totally asymmetric exclusion model, which 
has been solved exactly ||15|| . An atom in the immobile state changes to the running state 
if it succeeds in hopping to the right hand side, while a running atom becomes immobile if 
the site to the right is occupied by an atom in the immobile state. Such transitions between 
the immobile and running states of an atom model the effects of damping. LG Model A is 
thus characterized by the parameters p and a. LG Model B includes the possibility that the 
running atom at a randomly chosen site may change to the immobile state with probability 7 
prior to the motion of the atom takes place. This spontaneous convection from the running 
to the immobile state is supposed to be more important for weak driving forces. In order to 
capture the features of hysteresis, Braun et al. [|TJ|] went on to include in the basic model a 



majority rule in the conversion between the two states of the atoms. If the randomly chosen 
site is occupied by the leading atom in a compact group of r running atoms and it is blocked 
by a compact group of s immobile atoms in front, all the r + s atoms will turn into the 
immobile (running) state if r < s (r > s). These three models have the advantage of being 
easily implemented numerically by carrying out Monte Carlo simulations. 

It is useful to study similar models within the context of dynamical systems. A physical 
consideration is that parallel updatings of the states of the atoms may be more appropriate 



than the sequential updatings in the LG models studied in Ref. []14" |. In the present work, we 
propose analogous models with parallel updating rules in that all the atoms evolve in every 
time step according to updating rules, and hence the models become cellular automaton 
(CA). We performed numerical simulations on the models. Another advantage of casting 
the models in terms of parallel updating rules is that a more systematic analytical approach, 
analogous to those successfully applied to traffic flow models jl6lJl7l| ) may be applicable. Such 



approach focuses on the time evolution on the state of each of the sites, namely whether 
the site is occupied by an atom in the running or immobile state or unoccupied. In general, 
spatial correlations of gradually increasing spatial extent are introduced as time evolves. 
Equations can be written down relating the state of a site at time t + 1 to quantities at 
time t. By suitably decoupling the spatial correlations, a set of coupled nonlinear equations 
can be obtained with the fixed point corresponding to the solution in the long time limit. 
The complexity of the set of equations depends on the extent of spatial correlations retained 
after decoupling. The approach has the advantage that it gives the fraction of running atoms 
in the steady state together with other spatial correlation functions. As an illustration of 
the general idea of the approach, we study Models A and B with parallel updating rules 
and results are found to be in reasonable agreement with numerical simulations within 
the approximation of retaining correlations upto two sites. The present work, therefore, 



complements that of Braun et al. |T4| and suggests an alternative way of studying the 
various models proposed within the context of tribology. 

The plan of the paper is as follows. In Sec. II, we define the modified models with parallel 
updating rules and present the results obtained by numerical simulations. Section III reports 
results of our analytical calculations on model A and B. Results are compared with numerical 
simulations. We summarize the results in Sec. IV. 

II. MODELS AND NUMERICAL RESULTS 
A. Model A 

We consider N atoms on a one-dimensional lattice of L sites with periodic boundary 



condition. Following the basic lattice gas model |14[ , we modify the rules such that parallel 
updating is incorporated. The updating rules are: 

Rule Al: If the ith site on the lattice at time t is occupied by an immobile atom, it has a 
probability a to be in the advancing immobile state, i.e., state that flavours forward biasing 



5 



and a probability 1 — a to be in the retreating immobile state, i.e., state that flavours 
backward hopping . An advancing immobile atom can either hop into the empty (i + l)th 
site and become a running atom if the (i + 2)th site is not occupied by a retreating atom 
or be blocked by an atom occupying the (i + l)th site without changing its state, while a 
retreating atom can either hop into the empty (z — l)th site and stay in the immobile state if 
the (z — 2)th site is not occupied by a running or an advancing immobile atom or be blocked 
by an atom occupying the (z — l)th site without changing its state. 

Rule A2: If the zth site on the lattice at time t is occupied by a running atom, it can move 
to the empty (i + l)th site if the (i + 2)th site is not occupied by a retreating atom or stay 
at the zth site and remain in the running state if it is blocked by a running atom or change 
to the immobile state if it is blocked by an immobile atom at the (z + l)th site. 

Rule A3: If the zth site on the lattice at time t is empty and is sandwiched between a running 
or advancing immobile atom at the (z — l)th site and a retreating immobile atom at the 
(i + l)th site, then the atoms at the two neighboring sites are equally probable to hop into 
the zth site according to rules Al and A2. The atom that fails to hop at that time step will 
remain in its original state. 



Model A thus represents a modification of the LG Model A in Ref. []14j with parallel 
updating rules. The quantity of interest is the mobility B of the system which is defined 
as the fraction of atoms in the running state in the long time limit. We have carried out 
numerical simulations on Model A. Figure 1 shows the dependence of the mobility B as a 
function of the drift parameter a for different concentrations of atoms p. Note that only the 
range 1/2 < a < 1 corresponding to a biased field to the right hand side is shown, although 
the range < a < 1/2 can also be studied taking Model A as a CA model in its own right. 
To illustrate the basic features of the model, we have performed numerical simulations on 
systems with L = 1000. Typically, about 2000 time steps are sufficient for approaching 
the long time limit. The mobility B is obtained simply by counting the fraction of running 
atoms in the long time limit. For a < 1, an average over 50 random initial configurations 
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are performed. For a > 1/2, B = 1 for p < 0.32. For 0.32 < p < 1/2, the mobility B 
becomes unity at some critical value a c (p). For p > 1/2, the concentration is sufficiently 
high so that B < 1 for all values of a. 

The particular point of a = 1 deserves further discussion. It is found that for p > 1/2, 
the mobility B at a = 1 in the long time limit depends on the initial condition. The results 
at a = 1 shown in Fig.l correspond to B = (1 — p)/p, which are obtained by using the initial 
configuration in which all the atoms are immobile. For arbitrary initial configurations at 
a = 1, B is found to lie within the range (1 — p)/p < B <1. Similar results are obtained in 
our analytical approach, discussed in the next section, by treating the model as a dynamical 
system. 

B. Model B 

Model A forms the basic CA model for further modifications. In particular, while irre- 
versible transition into the running state for an isolated atom is strongly favourable in the 
high-field limit (a ~ 1), it is possible for a running atom to convert spontaneously to the 
immobile state in the weak field case. Following the LG model B in Ref. 0], we introduce 
the following rule in addition to the rules Al, A2, and A3 stated above: 

Rule Bl: Before the updating rules Al, A2, A3 are applied in each time step, every atom 
in the running state has a probability 7 to change its state to the immobile state and a 
probability 1 — 7 to remain in the running state. After this consideration, all the atoms on 
the lattice evolve according to the rules Al, A2, and A3. 

Rules Bl, Al, A2, and A3 define the CA Model B. While a tends to lead to a larger 
fraction of running atom, the parameter 7 counteracts the effect and tends to increase the 
number of atoms in the immobile state. Hence, the mobility B is generally lower for 7 7^ 
cases than the 7 = case for the same value of a. Figure 2 shows the values of B as a function 
of a for different values of 7 with the concentration fixed at p = 0.4, which corresponds to 
a concentration at which the atoms are isolated if they are uniformly distributed on the 
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lattice. Note that the mobilities converge to unity at a = 1 for different values of 7. Only 
when 7 = will the mobility becomes unity for a c < a < 1. The results shown are typical 
for p < 1/2. Figure 3 shows the results for p = 0.6, which corresponds to a concentration at 
which there are always some atoms with nearest neighbors if they are uniformly distributed 
on the lattice. In this case, B < 1 for all values of 1/2 < a < 1 and < 7 < 1. The 
mobilities converge to the same value for different values of 7 at a = 1. For 7 = 0, Model B 
reduces to Model A and the mobility B lies in the range (1 — p)/p < B < 1 with the precise 
value depending on the initial condition. 



C. Model C 



Following the LG model C in Ref. ||14|| , further modifications can be made by taking into 
account the influence of the state of the surrounding atoms on that of a single atom, i.e., the 
'crowding effect'. The modifications involve the considerations of the jamming of a running 
block of atoms (i.e., a compact group of nearest neighboring atoms in the running state) by 
an immobile block of atoms (i.e., a compact group of atoms in the immobile state). If these 
two adjacent blocks of atoms are sandwiched between two empty sites at the two ends, the 
state of all the atoms will then follow that of the larger block. This rule will be referred to 
as the majority rule. Furthermore, since the spontaneous transition of atoms in the running 
state to the immobile state should be suppressed by an increasingly stronger forward bias, 
the parameter 7 should be a-dependent. We impose the relation |L4| 7 = 70(1 — a) 2 on the 



parameter 7, where 70 is a model parameter. Hence, the updating rules for Model C can be 
state explicitly as below: 

Rule CI: At a certain time step, if the (i — r + l)th site to the ith site are all occupied by 
running atoms and the (i + l)th site to the (i + s)th site are all occupied by immobile atoms 
together with the condition that the (i — r)th site and the (i + s + l)th site are empty, then 
in the case of r > s (r < s), all the r + s atoms become running (immobile). Immediately 
after the changes, the states of the sites are then updated according to the rules of Model 
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B with 7 = 70(1 — a) 2 . 

Rule CI together with Model B define the CA Model C. Figure 4 shows typical results 
for the mobility B as a function of a for different values of p with 70 taking on a value 
close to 0. It is observed that the results look very similar to that of Model A except for 
the presence of hysteresis, a characteristic feature of the FK model [13], for p < 1/2. The 



values of B obtained by gradually increasing a from a = 0.5 to a = 1 are generally smaller 
than those obtained by gradually decreasing a from a = 1, and the difference between the 
mobilities for increasing and decreasing bias fields at a particular value of a increases with 



a as observed in the LG model in Ref. [Tj] 



III. ANALYTICAL APPROACH 

The CA models with parallel updating rules have the advantage that they can be treated 
analytically within the context of dynamical systems. The general idea is to establish the 
time evolution equations for the state on each site of the lattice. The equations, in gen- 
eral, involve spatial correlation functions. With suitable approximations typically involving 
proper decoupling of the correlations, a closed set of dynamical equations can be obtained. 
Such a set of equations can be treated as a dynamical mapping between quantities at time 
t + 1 and those at time t. Hence, following standard approaches in dynamical systems, the 
solution in the long time limit can be found by studying the fixed points (attractors) of the 
set of equations. Such an approach has been successfully developed for traffic flow models 
T6| , [T7H in which the cars, which is analogous to the atoms in the present models, can only 



move in one direction without backward diffusion movements. To illustrate the idea, we ap- 
ply the approach to study the modified CA models. It turns out that for Models A and B, 
reasonably good agreement with numerical simulations can be obtained by retaining spatial 
correlations involving two neighbouring sites only. Application of the method to Model C is 
difficult due to the built-in long spatial correlations in the majority rule of the model, and 
hence results are only reported for Models A and B. 
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A. Model A 



Although one can treat Model B directly and obtain results of Model A by setting 7 = 0, 
it is, however, illustrative to treat the simpler Model A first. We denote the states of the 
ith site at time t by the following set of boolean variables: 

1, if the i th site at time t is occupied by a running atom 

Ri(t) = { (i) 

0, otherwise 



h{t) 



1, if the i th site at time t is occupied by an immobile atom 

(2) 

0, otherwise 



1, if the i th site at time t is occupied 

(3) 

0, otherwise. 



Si(t) = Ri(t) + Ii(t) ~- 

Obviously, these variables satisfy the relationships Ri(t)Ri(t) = Ri(t)Si(t) = Ri(t), 



Ii(t)Ii(t) = Ii(t)Si(t) = hit), Si(t)Si(t) = 5,(0, Ri(t)Ri(t) = h(t)h(t) = 0, and 



Ri(t)Ii(t) = Ri(t)Si(t) = Ii(t)Si(t) = 0. Here Riit) represents the conjugate to Riit) given 



by Ri(t) — 1 — Ri(t), with similar definitions for Si(t) and Ii(t). 

In order to represent whether an immobile atom is advancing or retreating at a certain 
time step t, we define a boolean variable ijt (f) at the ith site at time t such that 



1, with probability / 
0, with probability (1 — /) 



M/) = < 

Thus the term 9 i>t (a)Ii(t) represents the probability that the ith site at time t is occupied 



by an advancing immobile atom, while the term 9 ijt (a)Ii(t) represents the probability that 



the ith site at time t is occupied by a retreating atom. Here 8^t(ct) denotes the conjugate of 



6> M (a), i.e., 9 ijt (a) = 1 -9 ijt (a). It follows that ijt (a)6 ijt (a) = 9 i:t (a), and 9 i:t (a)9 i:t (a) = 0. 
Similarly, in order to represent whether an advancing or an immobile atom can hop into 
the empty site sandwiched between the two atoms at a certain time step, another boolean 
variable r]i^(f) defined exactly the same as 9ij(f) is introduced with / = 1/2. With this, 



the factor r] itt (l/2)R i (t)S i+ i(t)[9i + 2 t t(oi)Ii+2(t)] represents, for example, the probability that 
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the advancing running atom at the ith site can hop successfully into the (i + l)th site at 
time t. Note that the two boolean variables 9 ijt (f) and rj ijt (f) are statistically uncorrelated. 

We study the time evolution of the variables Rt(t) and Ii(t), i.e., we seek the variables 
Ri(t + 1) and Ii(t + 1) as a function of quantities at time t. Focusing on Ri(t + 1), there are 
various ways in which the situation at time t affects Ri(t + 1). From the rules Al and A2, a 
running or an advancing immobile atom occupying the (i — l)th site at time t will hop into 
the empty ith site if the (i + l)th site is not occupied by a retreating immobile atom. At 
the next time step, the (i — l)th site will become empty while the ith site will be occupied 
by a running atom. This leads to a contribution to Ri(t + 1) of the form 

[Ri-i(t) +e i ^ t (a)I i _ 1 (t)][^][e i+1 , t (a)I i+1 (t) +S^j + Ri+^t)}. 

The three brackets express the conditions on the (i — l)th, ith, and (i + l)th sites at time 
t, respectively. The three terms in the last bracket is equivalent to saying that the (i + l)th 
site is not occupied by a retreating immobile atom at time t. 

The second contribution to Ri(t + 1) comes from the situation in which the (i + l)th 
site is occupied by a retreating immobile atom. In this case, the rule A3 leads to another 
probabilistic event. This situation contributes a term to Ri(t + 1) of the form 

^^[Ri^ii] + i _ 1)t (a)/ i _i(*)][Si(i)][0 i+ i it (a)/ i+ i(t)], 

which denotes the probability that the atom at the (i — l)th site succeeded in moving forward 
into the empty ith site and became a running atom. The first factor follows from 

the rule A3 as the atoms at the (i — l)th and (i + l)th sites are equally probable to hop into 
the ith site. 

Another contribution comes in when the ith site is occupied by a running atom, the 
(i + l)th site is empty and the (i + 2)th site is occupied by a retreating immobile atom and 
that the retreating atom succeeded in hoping back onto the (i + l)th site in the process. 
In this situation, the running atom will stay on the ith site at the next time step. This 
contributes a term 
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Vht (-)[R^[S l+1 (t)m + 2A^ + 2(t)] 

to Ri(t + 1), with the brackets expressing the conditions at the ith, (i + l)th, and (i + l)th 
sites. A fourth contribution comes from the situation that a running atom at the ith site is 
blocked by another running atom at the (i + l)th site and it gives a term 

Ri(t)R i+ i(t) 

to Ri(t + 1) according to the rule A2. 

Collecting all the four contributions to Ri(t + 1), we have 

Ri(t + 1) = [i2 + 0(a)/] i _ M [^| iit [0(a)/ + ^ + i2] i+M 

+lv(\)(R + e(a)i)] t ^ t [s^AW)iUi,t 
^v(\)R}iMi + iAW)i}i + 2,t + [R]i,t[R]i + i,u (5) 

which is a time evolution equation for Ri(t + 1) in that all the quantities on the right hand 
side are evaluated at time t. Note that we have simplified the notations so that all the 
quantities inside a squared bracket are to be evaluated at the position and time indicated 
as subscripts outside the bracket. 

Similar argument can be carried out for h(t + 1), although the analysis is slightly more 
complicated than the case of Ri(t + 1). The rule Al states that if the (i — l)th site is not 
occupied by a running or an advancing immobile atom, then the retreating immobile atom 
at the (i + l)th site can hop into the empty ith site deterministically. This contributes to 
Ii(t + 1) a term of the form 

[0~^/i-i(t) + S~Jt)] iW)] [e i+ i,t(a)Ii+i (t)} , 

where the terms in the first squared bracket is equivalent to saying that the (i — l)th site is 
not occupied by a running or an advancing immobile atom. 

The rule A3 comes into the consideration through various situations. If the (i — l)th site 
is occupied by a running or an advancing immobile atom, the retreating immobile atom at 
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the (i + l)th site still has half a chance to hop into the empty ith site, and contributes a 
term 



to Ii(t + 1). Note that when an advancing immobile atom fails to hop forward, it will stay 
at the site and remain immobile. In this case, an advancing immobile atom located at the 
ith site at time t will still be there at time t + 1 contributing a term to Ii(t + 1) of the form 



An analogous situation arises when the running or advancing immobile atom located at the 
(i — 2)th site succeeded in hopping into the empty (i — l)th site, leaving a retreating immobile 
atom at the ith site. This contributes a term to hit + 1) as 



Blocking by atoms in the nearest neighboring sites contributes the following terms. If an 
advancing immobile atom at the ith site is blocked by the an atom at the (i + l)th site, or 
a retreating immobile atom is blocked by an atom at the (i — l)th site, the state of the ith 
site at time t + 1 remains to be immobile. These two terms in 7j(i + 1) are represented by 
[6 i>t (a)Ii(t)]S i+1 (t) + Si_i(t)[Mo)Ji(*)]- Finally according to rule A2, a term Ri(t)I i+l (t) 
in Ii(t + 1) arises from the blocking of a running atom by an immobile atom. 

Collecting all the contributions to Ii(t + 1), we have 



77 i _i lt (-)[i2 i _i(*) + i _i lt (a)/ i _i(*)][5 , i (*)][0 i+ i )t (a)/ i+ i(t)] 



Vht {-)6 ht {a)I t {t) [S i+1 (t)][9 i+2tt (a)I i+2 (t)}. 



h(t + 1) = [0(a) J + S + V (-)(R + 0(a)/)] i _i it [^| M [0(a)7] i+lit 




+[^(i)(i2 + fl(a)i)] i _ 2>t [5| i _ 1>t [fl(a)i] i , t 



+[0(a)/] M [Sl i+ i )t + [Sl i _i, t [0(a)7] iit 



+[i*kt[JW 



(6) 
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Equations (§) and (||) can be used to compute the time evolution of the mobility of the system 
and the spatial averages of the products of different combinations of the state variables 
defined on the same or neighboring sites. 

The mobility B(t) at time t, i.e., the fraction of atoms in the running state at time t, 
can be expressed in terms of Ri(t) as 

B(t) = ±j:R>(t) = - o mt)}, (7) 

i i 

where (■■■) = l/NJ2i{~ ' ') is the spatial average of the quantity concerned over the system. 
It follows that B{t + 1) = l {Ri{t + 1)). Making use of the expression for Ri(t + 1) in terms 
of quantities at time t given by Eq.([|), the mobility at time t + 1 can be expressed in terms 
of spatial averages involving strings of upto three neighboring sites at time t. This gives 

B{t + 1) = [a(R0I) t + a 2 (I0I) t + (R00) t + a(I00) t + (R0R) t + a(I0R) t ] 

+±^[(R0I) t + a (I0I) t ] + L^(R0I) t + (RR) t . (8) 

For simplicity, we write the spatial averages at time t as (■ ■ -)t and express the strings of 
neighboring sites in order from left to right. We use the symbol "0" to denote an empty site 
or S. For example, {R0I) t implies counting the strings of neighboring sites with a running 
atom on the left and an immobile atom on the right with an empty site in between over 
the system at time t. The prefactors resulted from the fact that the spatial averages of the 
boolean variable 6{a) (or its conjugate) survives with a probability a (or 1 — a). Similarly, 
?7(l/2) survives with probability 1/2 under averaging. Thus, the term a(R0I) t comes from 



the average (RQ6(a)I) t . Noting that Si(t) + Si(t) = Ri(t) + kit) + S^t) = 1, we have 

(9) 



B(t + 1) = - 
P 



I- a 

(RR)t + (R0) t + a(I0) t - a — - — (IOI)i 



Equation (|9]) is an exact expression for B{t + 1) in terms of quantities evaluated at time 
t. In order to proceed, we write down the evolution equation for the spatial averages on 
the right hand side of Eq.(Q). Obviously iterating the equations backward in time gives 
terms involving longer strings of neighboring sites and hence longer spatial correlations. 
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To close the set of equations, a decoupling scheme retaining spatial averages involving two 
neighboring sites is invoked. The set of equations can then be treated as a dynamical system. 
The fixed points of the equations then give the results corresponding to the long time limit. 

To treat the system analytically, we decouple the term (IOI)t in Eq.(|9|) into products of 
averages involving two neighboring sites, i.e., (I0I) t ~ (I0) t {0I) t /(l — p), where 1 — p = (Q) t 
is the probability of finding an empty site |18| . With this approximation, Eq.(§) becomes 



B(t+1) = - 
P 



(RR^ + mt + ailOh-a 1 -^ 10 ^ 



(10) 



2 1 - p 

With the four variables R, I, S and S, a total of sixteen two-site spatial averages can 
be formed, among which four of them can be chosen to be independent. We choose the 
independent spatial averages to be (RR)t, (RI)t, (IR)t an d (H)t- The other two-site averages 
are related through 

(Rl) t = (RR) t + (RI) t , 
(lR) t = (RR) t + (IR) t , 
(Il) t = (IR) t + (II) t , 
(U) t = (RI) t + (J/) t , 
(R0) t = P B{t) - (Rl) t , 
(I0) t = p(l-B(t))-(Il) t , 
(0R)t = P B{t) - (lR) t , 
(0I) t = p(l-B(t))-(lI) t , 
(01) t = p-(Rl) t -(Il) t , 

(10) t = (01) t , 

(00) t = l-p-(01) t , 

(11) t =(/l) t + (i21) t , (11) 

where "1" represents an occupied site regardless the character of the atom. Using (RR) t+1 = 
(Ri(t + l)Ri +1 (t + 1)) and Eq.(|) for R^t + 1) and R i+1 (t + 1), we have 
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{RR) t+1 = (RRR) t + 



1 — a 



(RROI)t + (RORR) t + a(IORR) t 



(12) 



To make the approximation self-consistent, we invoke the decoupling scheme and retaining 
spatial averages involving no more than two sites, we have 



(RR) 



t+i 



+ 



(RR)t 



(RR) 2 t 

pB(t) 
l-a(R0) t {0R) t (0I) t 



1 — a 



l-p) P B(t) I 2 

S '[(R0) t + a(I0) t 



(R0){0I) t +{0R) t [{R0) t + a{I0) t ]} 



4 2 (l-pf P B(t) 
Similarly, for (RI) t+1 , (IR) t+1 , and (77) t+1 , we obtain after decoupling 

1 



(RI) 



t+i 



(RR)t{Rl)t | 1 - a {QI)t 



a 



(RD) t + -(I0) t 



(00) t (0I) t + a 



(IR) 



t+i 



pB(t) 1-p 
j (0R) r (RI) t I- a 
\ pB(t) + 1-p 

a l-a (0I) 2 t (I0) t \ 
+ 21-pp(l-B(t))J' 
1 - a 

(RR) t + — -(R0) t (I0) t 



+ 



1-p 
p(l-B(t)) 



[{RQ) t + a(I0) t ] 



pB(t) 



2(1 -P) 



a(IR) t + (l-a) 



(u)t(iR)t 

p(l-B(t)) 



a (0I) t (IR)t 



2(1 - p) p(l - B(t)) 



(II) t+1 = (1 - a)[(RI) t + a(II) t ] + ^L-^(0/> t (/0> t + ^-ARI) t (IR) t 



2 (lI) t (II) t 
+ {1 ~ a) P (l-B(t)) +a 
| (1 - af (01} t (ll) t (10} t | 



2 1 - p 



a- 



P B(t) 



p(l-B(t)) 1 
l-a(RI) t (U)t(IR)t 



1-p p(l-B(t)) pB(t) p(l-B(t)) 

4W1 (/0> t (0i) t 
+a(l - aj^TT^ ^ /jXN o + 



.[(i2/) t + a(//) t ] + 



p 2 (l-B(t)f 2 l-pp 2 (l-B(t)) 2 
1-a) 2 (0I) t (II) t 



2(1 -£}(*)) 



[(R0)t + a(I0) t 



I I- a 



(07)1(70), 



2 V 1 -p/ p(i -£(*)) 

(07) t (17) t (77) t (70) f | 1-a 



[(R0) t + a(I0) t ] + 



a (1 — a)' 



2 1-p 
(0I)t(IR) t (RI) 



2(1 -p) P B(t) p(l-B(t)) 



+a{I0) t ] + 



al-a(71) t (07) t (77), 
2 1-p p 2 (l -B(t)) 2 



[{R0) t + a{I0) t ] 



(13) 



(14) 



(15) 
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+ 4iwj p\l-B{t)f [W t + «(/0) t ]. (16) 

Equations flUf ) and (|i~3|)- ( pT6|) form a set of five equations for £>, (RR), (RI), (IR), 
and (//)• These equations form a five dimensional dynamical system. To compare with 
simulation data, we solve for the stable fixed points numerically. Results for the mobility 
are shown in Fig.l as solid lines for different values of a and p. The analytical results are in 
good agreement with numerical data, showing that the decoupling approximation is sufficient 
to capture the essential features of the model. We have also checked the results of the two- 
site and three-site spatial averages numerically and analytically and it is found that the 
decoupling scheme gives qualitatively correct results for the spatial averages. Our method 
represents a systematic way of deriving mean field theories from microscopic consideration 
by following the time evolution of the system. The decoupling scheme of retaining two-site 
spatial averages is the minimal procedure to obtain qualitatively correct mobility and spatial 



averages involving longer strings of sites fi~8|j 



The a = 1 case deserves further discussion. For a = 1, Eq.(|9]) gives 

B(t + l) = - p [(10) t + (RR) t ], (17) 

where (10) t = (R0) t + (I0)t- The spatial average (RR) t +i can be obtained by setting a = 1 
in Eq.([L3|) to get 

(10) t (0R) t (RR) t (RR)* f . 

{RR)t+1= a-pwt) (18) 

To form a closed set of equations, we work out the spatial averages (10) t+ i and (0R) t +i 
within the approximation of retaining two-site correlations to get 

<10> t = p-(l-^)(p-<10> t + Mi), (19) 

p l-p 



and 



/nm /inU !Mi (W) t (0R) t (RR)t (9n , 
(0R) t+1 = (10), + — — { i-p)p B (t) ' (20) 
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where we have used the relations stated in Eg. (|TT|) . The fixed points satisfy B{t + 1) = 
B(t) = B, (10)*+! = (10) t = y, (RR) t+ i = (RR) t = z, and (0R) t+1 = (0R) t = w. Hence, 
they satisfy the simltaneous equations 

P B = y + z, (21) 

V = P-(1--)(P-V + T?-), (22) 

z 2 wyz .„„. 

wyz . jS 

w = y + ^ . (24) 

p£ (l-p)p5 1 J 

From Eg. (p2|) , the stable fixed point for y is given by y = p for p < 1/2 and y = 1 — p for 
p > 1/2. From Eq. (|23|) , z = is a stable fixed point for p < 1/2 and z = pB — w is a stable 
fixed point for p > 1/2, where we have used the results for y. It is important to note that 
Eq.(p4|) for w becomes redundant. Hence the situation is that we have three equations with 
four unknowns. The values of B and z are governed by the linear relation 

z = pB-(l-p). (25) 

Any values of x G [(1 — p)/p, 1] and z G [0,2p — 1] satisfying Eq.(p5[) is a solution to the 
system of equations. Thus for a = 1, the mobility 5 = lforp<l/2 and B lies in the 
range [(1 — p)/p, 1] for p > 1/2 with the precise value depending on the initial condition, in 
agreement with numerical results. The value (1 — p)/p shown in Fig.l corresponds to the 
initial condition of all the atoms being immobile. 



B. Model B 

The new parameter 7 introduced in rule Bl is the probability that an atom in the running 
state changes into the immobile state in a time step. To carry out analytical treatments 
similar to those in Model A, it is convenient to divide each time interval into two halves. In 
the first half of a time step, rule Bl applies and the parameter 7 is effective; while in the 
second half of the time step, rules Al and A2 apply. Introducing a Boolean variable Cm(t) 
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analogous to, but statistically independent of, 9i,t{l) and rj^t, the variables Ri(t) and Ii(t) 
evolve in the first half of the time step as: 

Ht + l) = QMRi(t), (26) 

and 

i i {t + ±) = i i {t) + &MRi(t)- (27) 

The time evolution in the second half of the time step is given by Eqs. @ and @, with 
the quantities on the right hand side of the equations corresponding to those evaluated at 
t + 1/2. Combining the evolution in the two halves of a time step, we finally arrive at 

Riit + 1) = {[CvT) + *(a)C(7)]fl + ^a)/}*-!,^]^ + (C[t) 

+fl(a)C(7))i2 + »(a) J] m ,t + {v(\)iOl) + 0^)C(l)]R 

+ [C^)R] t , t [W)RWi,t, (28) 

and 

hit + 1) = (S + [C(7)^R + rK~)C(7) + »7(^(a)C(7)] J? 
+p^ + i7(i)e(a)]J}^x,t[5| 4 , t pRC(7)i2 
+9{a)I} i+1>t + [^(i)e(a)C(7)i2 + ?K~) Wk* 

^ i+ iAW^ai)R+W)ni+2, t +{[v(l)aj) 

+v(\)d(<x)((l)]R + v(l)0(a)Ih- 2 ASU,t 

i9^C(i)R + W)ih + [d(a)C(i)R 

+e(a)I] i , t [S} i+1 , t + [S]^ ltt [e{a)C(l)R 

+9{a)Ih + [W)Rh[((l)R + W (29) 
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Equations ( p8|) and (p9|) are the time evolution equations relating + 1) and Jj(t + 1) to 
quantities at time t. They play exactly the same role as Eqs. (^|) and (|6|) in model A. 

It is then straightforward to carry out the same treatment for model B as in model A, 
and we simply outline the key steps in the following discussion. Following the same steps 
leading to Eq. ([]), the mobility B(t + 1) at time t + 1 for model B is given by 

B(t + 1) = ~{(1 - i) 2 (RR) t + (1 - 7 + ory){R% + a(I0) t 
-|(1 - a)[(I0I) t + i((R0I) t + (I0R) t ) 
W(R0R)t]}- (30) 

Equation (|30"D is the generalization of Eq.@ to model B. It reduces to Eq.@ for 7 = 0. 
Employing a decoupling approximation to retain only spatial averages involving upto two 
nearest neighboring sites, Eq. (|30|) becomes 

B{t + 1) = -{(1 - !) 2 {RR) t + (1 - 7 + <x-y)(R0) t + a(I0) t 
P 

"f T^7 [(/0)t + ^ R0)t ^ 0I)t + 7(0jR) *] } - (31) 

To close the set of equations, we construct the time evolution equations for the spatial 
averages {RR), (RI), (IR), an d (//). The other spatial averages can be constructed from 
these four averages. In the presence of the parameter 7, the resultant equations are more 
complicated than Eqs.(13)-(16) in model A. This set of equations forms a dynamical system. 
The stable fixed point corresponding to the mobility and spatial averages in the long time 
limit can be readily solved numerically. Results for the mobility B in the steady state 
are shown as solid lines in Figs. 2 and 3 for two different values of atomic concentration 
p. The theoretical results obtained within the decoupling approximation capture all the 
essential features of the numerical data. It is observed that the theoretical results are 
consistently slightly greater than the numerical data. The discrepancies come from the 
decoupling scheme. If the decoupling approximation is extended to retain spatial averages 
involving upto three neighboring sites, for which the calculations are much more involved, 
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the results are in better agreement with numerical data |19[ . It is, however, important to 
stress that the essential physics is captured within the two-site decoupling approximation. 

The particular case of a = 1 can be treated in a way analogous to that in model A. It is 
found that for 7 e (0, 1], the stable attractors give the mobility 

I 1 if p < 1/2 

B{t - 00) = , (32) 

( ^ if p > 1/2 

together with the spatial averages 



(10) t -oo = (0i?)^oo 

p if p < 1/2 
1 - p if p > 1/2 



(33) 



and 



(RR)^ = 0. (34) 

For a = 1 and 7 7^ 0, the time evolution of the state at a site depends on the states of 
the nearest neighboring sites only and hence the decoupling scheme retaining only two-site 
spatial averages is good. Results so obtained are in exact agreement with numerical data. 
It should be noted that for 7 = and a — 1, B — 1 for p < 1/2 and B lies in the range 
[(1— p)/p, 1] for p > 1/2 with the precise value depending on the initial condition as discussed 
in the previous subsection. 



IV. SUMMARY 

In summary, we have proposed three CA models defined in terms of parallel updating 
rules analogous to the three models recently studied by Braun et al. which are defined 
in terms of sequentially updating rules. These models are of increasing complexity so as to 
model the generalized FK models proposed recently within the context of tribology. The 
first model (Model A) involves atoms in two different dynamical states, i.e., running and 
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immobile, subjected to an external field parametrized by a. Atoms in the running state tend 
to hop along the field direction while atoms in the immobile state may bounce backward. 
The second model (Model B) involves spontaneous transition of atoms from the running state 
to the immobile state in addition to the rules of Model A. The third model (Model C) takes 
into account of the crowding effect of the system as well. Results of numerical simulations 
indicate that the mobility, which is defined as the fraction of atoms in the running state, as 
a function of a for the three models exhibit phase transitions, jammings and hysteresis. The 
proposed CA models have the advantage that they can be treated as dynamical systems and 
analysed by the microscopic approach previously applied to analogous models of traffic flow 
problems ||17|| . For Models A and B in which the state of a site is affected by the nearest and 
next nearest neighbours, time evolution equations of the state of the sites can be explicitly 
written down. The mobility of the system at time t + 1 can be expressed in terms of spatial 
averages at time t involving not more than three sites. By invoking an approximation 
which decouples the spatial averages into products of spatial averages involving upto two 
neighboring sites as well as deriving the evolution equations of the two-site spatial averages, 
a closed set of equations forming a nonlinear dynamical mapping can be established. The 
behavior in the long time limit can be found by studying the stable fixed point of the 
mapping. The fixed point of the mapping can be found numerically. In the special case that 
a = 1 and 7 ^ 0, the decoupling scheme is exact and analytical solutions can be found. For 
the whole range of possible values of the parameters in the models, the present approach 
yields satisfactory results when compared with simulation data. Our analytical approach 
represents a systematic way to derive approximations of increasing accuracy starting from 
a microscopic point of view by following the time evolution of the system. It is expected 
that better agreement with numerical results can be obtained by retaining spatial averages 
involving a long string of neighboring sites. 
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FIGURES 

FIG. 1. The mobility B, which is defined as the fraction of atoms in the running state, of Model 
A in the long time limit as a function of the external biasing parameter a for different values of the 
concentration of atoms p. The symbols are numerical data and the solid lines are results obtained 
by invoking the decoupling scheme in the analytical approach discussed in the text. 

FIG. 2. The mobility B of Model B in the long time limit as a function of the parameter a 
for different values of the parameter 7 at fixed concentration p = 0.4. The symbols are numerical 
results and the solid lines are results obtained by invoking the decoupling scheme. Results are 
typical of those for p < 1/2. 

FIG. 3. The mobility B of Model B in the long time limit as a function of the parameter a 
for different values of the parameter 7 at fixed concentration p = 0.6. The symbols are numerical 
results and the solid lines are results obtained by invoking the decoupling scheme. Results are 
typical of those for p > 1/2. 

FIG. 4. The mobility B of Model C in the long time limit as a function of the parameter a for 
different values of the concentration p. The parameter 70 is taken to be 10~ 5 . Results are typical 
for those with 70 ~ 0. Hysteresis in the mobility is observed for p < 1/2. 
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